## --------------------------------------- ##
## simulation setup:
##  delta = n / p
##  rho = k / n
## --------------------------------------- ##

## --------------------------------------- ##
## cluster setup
## --------------------------------------- ##
## i <- as.numeric(Sys.getenv("SLURM_ARRAY_TASK_ID")) # this should be 1--50
for(i in 1:length.out){
  cl <- makeCluster(detectCores()/2)
  registerDoParallel(cl)
  

    ## draw error term and n
    set.seed(38193)
    n    <- floor(p * delta[i])
    eta  <- rnorm(n)
    etaz <- sigma * eta
    
    ## doRNGseed(1234)
    sim_out <- foreach(j = 1:length(rho), .combine = "rbind", .packages= c("BridgeChange", "glmnet")) %dorng% ({
        
        cat("Now at i = ", i, " j = ", j ,"\n")
        k <- ceiling(n * rho[j])
        
        l2_error <- save_bridge <- save_ridge <- save_en  <- rep(NA, n_iter)
        for (iter in 1:n_iter) {
            rm(y); rm(X);
            ## gen data ***************************************************************
            X <<- matrix(rnorm(n * p), nrow = n, ncol = p)
            beta <- c(rnorm(k, 0, 3), rep(0, p - k))
            y <<- X %*% beta + etaz
            cat("dim(X) = ", dim(X), " length of y ", length(y), "\n")
            ## fit lasso ***************************************************************
            ## out <- rlasso(y = y, x = X, intercept = FALSE)
            ## l2_error[iter] <- sqrt(sum((coef(out) - beta)^2)) / sqrt(sum(beta^2))
            out <- cv.glmnet(y = y, x= X, alpha = 1, nfolds = 3)
            l2_error[iter] <- mean((X %*% (coef(out, s = "lambda.min")[-1,] - beta))^2)
            
            ## fit BridgeChange ********************************************************
            ## out1 <- BridgeChangeReg(y=y, X= X, scale.data=TRUE, intercept=TRUE,
            out1 <- BridgeChangeReg(y~X, mcmc= 100, burn = 100, thin=1, verbose=0, n.break = 0)
            beta.bridge <- coef_bridge(out1)
            save_bridge[iter]  <- mean((X %*% (beta - beta.bridge))^2)
            
            ## fit ridge ********************************************************
            rg_out <- cv.glmnet(y = y, x= X, alpha = 0, nfolds = 3)
            save_ridge[iter] <- mean((X %*% (coef(rg_out, s = "lambda.min")[-1,] - beta))^2)
            
            ## fit elastic net **************************************************
            en_out <- cv.glmnet(y = y, x= X, alpha = 0.5, nfolds = 3)
            save_en[iter] <- mean((X %*% (coef(en_out, s = "lambda.min")[-1,] - beta))^2)
            
            
        }
        
        mse.lasso   <- median(l2_error, na.rm = TRUE)
        mse.bridge  <- median(save_bridge, na.rm = TRUE)
        mse.ridge   <- median(save_ridge, na.rm = TRUE)
        mse.elastic <- median(save_en, na.rm = TRUE)
        
        ## output
        out <- c(mse.lasso, mse.bridge, mse.ridge, mse.elastic)
        out
    })
    cat("\nCurrent simulation is ", i, " iteration.\n")
    saveRDS(sim_out, file = paste("./nochange/corr00/res_pred/sim_pred_mse_corr00_", i, ".rds", sep = ''))
    
    
    stopCluster(cl)
    stopImplicitCluster()
}
